Tarea 4: Bayesian Inference

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

Bienvenides a la cuarta tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

#Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación y funcional
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tidyr)
library(purrr)

# Plots
library(scatterplot3d)
## Warning: package 'scatterplot3d' was built under R version 4.5.2
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
# Varios plots
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
# Análisis bayesiano
library(rethinking)
## Loading required package: rstan
## Warning: package 'rstan' was built under R version 4.5.2
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 4.5.2
## 
## rstan version 2.32.7 (Stan version 2.32.2)
## 
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## 
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## 
## Attaching package: 'rstan'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## Loading required package: parallel
## Loading required package: dagitty
## rethinking (Version 2.01)
## 
## Attaching package: 'rethinking'
## 
## The following object is masked from 'package:purrr':
## 
##     map
## 
## The following object is masked from 'package:stats':
## 
##     rstudent

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Por último, si sigue con problemas, intentar con el siguiente script:

https://github.com/dccuchile/CC6104/blob/master/scripts/rethinking_install.R

En caso de no lograr instalar la librería a pesar de estos pasos, favor contactar al auxiliar por correo.

Pregunta 1: Tasa de conversión en e-commerce (Binomial)

Se dispone del archivo local online_shoppers_intention-2.csv, donde la variable Revenue indica si un usuario realizó una compra (TRUE/FALSE o 1/0). Se busca estimar la tasa de conversión promedio y analizar cómo cambia la inferencia bayesiana al variar el prior y el tamaño de muestra.

Se pide:

  1. Cargue el archivo y convierta Revenue a variable binaria {0,1}. Calcule la tasa global de conversión.
    • ¿Qué proporción de usuarios realizó una compra?
  2. Implemente una grid approximation para \(p \in [0,1]\) con 1000 puntos usando los priors:
    • Prior A: \(p \sim N(0.14, 0.05)\)
    • Prior B: \(p \sim \text{Uniforme}(0,1)\).
  3. Para tamaños de muestra \(n \in \{10, 50, 500, 1000\}\), muestree datos aleatorios, calcule \(x=\sum \text{Revenue}\) y grafique las posteriores \(p\mid x,n\) para ambos priors y comparelo con el valor real.
    • ¿Cómo cambia la forma de la posterior al aumentar \(n\)?
  4. Ajuste el modelo con Laplace approximation (quap) bajo prior uniforme y compare sus resultados con el método de grilla.
    • ¿Los parámetros estimados son coherentes con los obtenidos por grid approximation?
  5. Calcule los intervalos de credibilidad y HPDI al 50% y 95%.

    • ¿Son similares estos intervalos? ¿Porqué?

    • ¿Que diferencia tiene con respecto al intervalo de confianza de el enfoque frecuentista?

Respuesta:

  1. Cargue el archivo y convierta Revenue a variable binaria {0,1}. Calcule la tasa global de conversión.
    • ¿Qué proporción de usuarios realizó una compra?
#Fijamos la semilla para mantener los mismos resultados
set.seed(3)

datos_ecom <- read_csv("datos/online_shoppers_intention-2.csv", show_col_types = FALSE)


datos_ecom$Revenue <- ifelse(datos_ecom$Revenue == "FALSE", 0, 1)

suma <- sum(datos_ecom$Revenue)
total <- nrow(datos_ecom)

prop <- suma/total

#Así está bien o en porcentaje?
cat("Con un total de ", total, " datos, y una cantidad ", suma, " de compras, ")
## Con un total de  12330  datos, y una cantidad  1908  de compras,
cat("la proporción de usuarios que sí realizó la compra es de: ", round(prop*100, 2), "%")
## la proporción de usuarios que sí realizó la compra es de:  15.47 %
  1. Implemente una grid approximation para \(p \in [0,1]\) con 1000 puntos usando los priors:
  1. Para tamaños de muestra \(n \in \{10, 50, 500, 1000\}\), muestree datos aleatorios, calcule \(x=\sum \text{Revenue}\) y grafique las posteriores \(p\mid x,n\) para ambos priors y comparelo con el valor real.

¿Cómo cambia la forma de la posterior al aumentar \(n\)?

#Los tamaños de muestras que nos piden
muestras <- c(10, 50, 500, 1000)


posterior <- function(prior, sec, nombre){
  par(mfrow=c(2,2))
  for (n in sec){
    muestra <- sample(datos_ecom$Revenue, size = n, replace = FALSE)
    x <- sum(muestra)
    
    #probabilidad de observar exactamente x éxitos en n muestras para cada una de las p probabilidades de p_grid
    likelihood <- dbinom(x, n, prob=p_grid)
    
    #posterior sin estandarizar
    unstd.posterior <- likelihood * prior
    
    #la estandarizamos y tenemos la distribución actualizada de p dada nuestros datos
    posteriors <- unstd.posterior / sum(unstd.posterior)
    
    max_post <- max(posteriors)
    max_index <- which.max(posteriors)
    max_p <- p_grid[max_index]  
    
    print(paste("Máximo de posterior para ", nombre, " con n = ", n, ": ", round(max_post, 4), " en p = ", round(max_p, 4)))
    
    plot(p_grid, posteriors, 
     xlab="Probabilidad de realizar la compra",
     ylab="Probabilidad posterior",
     main = paste("Posterior con prior", nombre , "- n =", n),
     xlim = c(0, 1), ylim = c(0, max(posteriors)))
    abline(v=prop, col="red", lwd=2, lty=2)
    
    legend("topright", legend=c("Posterior", "Proporción"),
         col=c("black","red"), lty=c(1,2), lwd=2)

    
  }
}

posterior(prior_A, muestras, "A")
## [1] "Máximo de posterior para  A  con n =  10 :  0.0088  en p =  0.1331"
## [1] "Máximo de posterior para  A  con n =  50 :  0.0115  en p =  0.1502"
## [1] "Máximo de posterior para  A  con n =  500 :  0.0265  en p =  0.1471"
## [1] "Máximo de posterior para  A  con n =  1000 :  0.0346  en p =  0.1712"

posterior(prior_B, muestras, "B")
## [1] "Máximo de posterior para  B  con n =  10 :  0.0029  en p =  0.3003"
## [1] "Máximo de posterior para  B  con n =  50 :  0.0104  en p =  0.0801"
## [1] "Máximo de posterior para  B  con n =  500 :  0.0246  en p =  0.1562"
## [1] "Máximo de posterior para  B  con n =  1000 :  0.0348  en p =  0.1562"

¿Cómo cambia la forma de la posterior al aumentar \(n\)?

Al ir aumentando n, la posterior se va pareciendo cada vez más a la posterior real. Esto se puede observar al ver que el máximo del gráfico se va acercando cada vez más a 0.15, la proporción que sacamos anteriormente. Esto se puede explicar porque a medida que aumenta la cantidad de datos muestreados, los datos predominan por sobre los priors, concentrando la curva alrededor de la proporción real.

Compararando ambos priors, el prior A varía menos a medida que aumenta n, convergiendo suavemente a la distribución real. Esto se puede explicar, dado que parte con una media en 0.14 muy cercana a la proporción real de unos 0.15. Mientras que el prior B, al inicio tiene menos certeza (todos los valores son igual de probables), entonces con menos datos tiene menos información para acercarse al valor real, pero a medida que aumenta n, va pareciéndose más a dicha proporción.

  1. Ajuste el modelo con Laplace approximation (quap) bajo prior uniforme y compare sus resultados con el método de grilla.
#Para poder comparar los cuatro gráficos de una
par(mfrow=c(2,2))

for (n in muestras){
  muestra <- sample(datos_ecom$Revenue, size = n, replace = TRUE)
  
  x <- sum(muestra)
  
  #Como lo vimos en clase
  globe.qa <- quap(
   alist(
   x ~ dbinom(n, p), # binomial likelihood
   #Por problemas con el knit cambiamos el dunif(0,1) por dbeta(2,2), que es lo más cercano a uniforme posible que no tenga problemas con x = 0 o x = n
   p ~ dbeta(2,2) # uniform prior
    ) ,
    data=list(x=x, n=n))
   
  # Extraemos el valor de p desde el globe.qa
  posterior_quap <- extract.samples(globe.qa)
  
  #Densidad de p obtenida de la aproximación
  densidad<- density(posterior_quap$p)  
 
  max_densidad <- max(densidad$y)
  max_p <- densidad$x[which.max(densidad$y)]
  print(paste("Máxima densidad de Laplace: ", round(max_densidad, 4), " en p = ", round(max_p, 4)))
  
  plot(densidad, main = paste("Posterior con Laplace - n =", n),
       xlab = "Probabilidad de compra", ylab = "Densidad",
       xlim = c(0, 1), ylim = c(0, max(densidad$y)))
      abline(v=prop, col="red", lwd=2, lty=2)
    
    legend("topright", legend=c("Posterior", "Proporción"),
         col=c("black","red"), lty=c(1,2), lwd=2)
  
}
## [1] "Máxima densidad de Laplace:  5.1134  en p =  0.0762"
## [1] "Máxima densidad de Laplace:  10.812  en p =  0.0832"
## [1] "Máxima densidad de Laplace:  25.0051  en p =  0.1424"
## [1] "Máxima densidad de Laplace:  34.8464  en p =  0.1527"

Al comparar Laplace con la grid approximation, ambos métodos producen resultados coherentes en general. Se observan algunas diferencias cuando n es pequeño y podría deberse principalmente a la variabilidad de la muestra, porque estamos analizando muestras aleatorias distintas. A pesar de esto, con n grandes (500 y 1000) la posterior se concentra y la aproximación de Laplace coincide muy bien con la solución por grid.

Veamos cómo se vería si comparamos con muestras iguales

comparacion <- function(){
  
  for (n in muestras){
    
    cat("Para n: ", n)
    
    # Usamos la misma muestra
    muestra <- sample(datos_ecom$Revenue, size = n, replace = FALSE)
    x <- sum(muestra)
    
    # grid
    likelihood <- dbinom(x, n, prob = p_grid)
    unstd_post <- likelihood * prior_B
    posteriors <- unstd_post / sum(unstd_post)
    
    # MAP grid
    MAP_grid <- p_grid[which.max(posteriors)]
    
    cat("MAP (Grid):    ", round(MAP_grid, 4), "\n")
    
    
    #laplace
    modelo <- quap(
      alist(
        x ~ dbinom(n, p),
        p ~ dunif(0,1)
      ),
      data = list(x = x, n = n)
    )
    
    muestras_quap <- extract.samples(modelo)
    dens_quap <- density(muestras_quap$p)
    
    # MAP
    MAP_laplace <- dens_quap$x[which.max(dens_quap$y)]
    
    cat("MAP (Laplace): ", round(MAP_laplace, 4), "\n")
  
    plot(
      p_grid, posteriors, type = "l", col = "blue",
      xlab = "p", ylab = "Posterior",
      main = paste("Grid vs Laplace (n =", n, ")")
    )
    lines(dens_quap$x, dens_quap$y / max(dens_quap$y) * max(posteriors), 
          col = "red", lwd = 2)
    
    legend("topright",
           legend = c("Grid", "Laplace"),
           col = c("blue", "red"), lwd = 2)
  }
}

comparacion()
## Para n:  10MAP (Grid):     0.3003 
## MAP (Laplace):  0.3044

## Para n:  50MAP (Grid):     0.1201 
## MAP (Laplace):  0.1206

## Para n:  500MAP (Grid):     0.1481 
## MAP (Laplace):  0.1514

## Para n:  1000MAP (Grid):     0.1572 
## MAP (Laplace):  0.158


Podemos observar que ahora ambas técnicas producen resultados aún más coherentes, donde los MAP son parecidos entre ambas y se parecen aún más cuando n aumenta.

  1. Calcule los intervalos de credibilidad y HPDI al 50% y 95%.
par(mfrow=c(2,2)) 

for (n in muestras){
  muestra <- sample(datos_ecom$Revenue, size = n, replace = TRUE)
  
  x <- sum(muestra)
  
  globe.qa <- quap(
    alist(
      x ~ dbinom(n, p),
      p ~ dunif(0, 1)   
    ),
    data=list(x=x, n=n)
  )
  
  posterior_quap <- extract.samples(globe.qa)
  
  cred_interval_50 <- PI(posterior_quap$p, prob=0.5)
  cred_interval_95 <- PI(posterior_quap$p, prob=0.95)
  
  hpdi_50 <- HPDI(posterior_quap$p, prob=0.5)
  hpdi_95 <- HPDI(posterior_quap$p, prob=0.95)
  
  cat("Intervalo de Credibilidad 50%: ", cred_interval_50, "\n")
  cat("Intervalo de Credibilidad 95%: ", cred_interval_95, "\n")
  cat("HPDI 50%: ", hpdi_50, "\n")
  cat("HPDI 95%: ", hpdi_95, "\n")
  
  densidad <- density(posterior_quap$p)
  plot(densidad, main = paste("Posterior con Laplace - n =", n),
       xlab = "Probabilidad de compra", ylab = "Densidad",
       xlim = c(0, 1), ylim = c(0, max(densidad$y))) 
  
  abline(v = cred_interval_50, col = "red", lwd = 2, lty = 2) 
  abline(v = cred_interval_95, col = "blue", lwd = 2, lty = 2) 
  abline(v = hpdi_50, col = "green", lwd = 2)  # HPDI 50%
  abline(v = hpdi_95, col = "purple", lwd = 2)  # HPDI 95%
  
    legend("topright",
         legend = c("PI 50%", "PI 95%", "HPDI 50%", "HPDI 95%"),
         col = c("red", "blue", "green", "purple"),
         lty = c(2, 2, 1, 1),
         lwd = 2,
         bty = "n",
         cex=0.7)
}
## Intervalo de Credibilidad 50%:  0.1166846 0.284856 
## Intervalo de Credibilidad 95%:  -0.04031261 0.4502999 
## HPDI 50%:  0.1159058 0.2837031 
## HPDI 95%:  -0.04045292 0.4493685
## Intervalo de Credibilidad 50%:  0.08907396 0.1514695 
## Intervalo de Credibilidad 95%:  0.0306897 0.2108611 
## HPDI 50%:  0.09439276 0.1562724 
## HPDI 95%:  0.02904603 0.2086294
## Intervalo de Credibilidad 50%:  0.1376426 0.1588014 
## Intervalo de Credibilidad 95%:  0.1165092 0.1788537 
## HPDI 50%:  0.1376402 0.1587902 
## HPDI 95%:  0.1177929 0.1798519
## Intervalo de Credibilidad 50%:  0.1444415 0.1596034 
## Intervalo de Credibilidad 95%:  0.1298629 0.1746387 
## HPDI 50%:  0.1438393 0.1589521 
## HPDI 95%:  0.1299458 0.1747096

A work by CC6104